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Abstract 

An  investigation  has  been  made  of  the  electron-electron  two  stream  instability 
as  the  beam  temperatures  were  increased,  using  particle  simulations.  Growth  rates 
and  saturation  characteristics  were  studied  and  compared  to  theoretical  models. 
The  final  state  of  the  system,  characterised  by  the  value  of  the  distribution  function 
at  zero  velocity,  /(v  =  0),  was  found  to  evolve  to  a  Maxwellian  in  the  cold  case,  and 
to  a  double-peaked  distribution  stable  by  the  Penrose  criterion  in  the  warm  case. 
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Amy  Wendt  (Prof.  C.  K.  Birdsall) 
May  3,  1985 


Introduction 


Rainer,  et.  al.  (1]  observed  a  transition  in  the  character  of  the  beam  plasma  in¬ 
stability  as  the  density  of  the  beam  was  varied  relative  to  that  of  the  plasma.  Our 
study  examines  the  changes  in  the  electron-electron  two  stream  instability  as  the 
beam  temperatures  are  increased,  in  search  of  a  similar  transition.  ESI,  a  particle 
simulation  computer  code,  was  run  with  beam  temperatures  varying  from  zero  to 
so  warm  that  no  instability  occurred.  In  particular,  the  program  was  used  to  study 
the  growth  rates  of  unstable  modes  and  observe  the  time  evolution  of  the  velocity 
distribution.  The  simulation  results  were  compared  with  theoretical  growth  rates 
and  stability  condiions.  No  well-defined  transition  as  observed  for  the  beam- plasma 
instability  was  found. 


Theory 


The  theoretical  growth  rates  of  the  various  modes  of  the  two  stream  instability  are 
given  by  the  dispersion  relation  D(u,  k)  —  0.  This  relation  can  be  found  by  solving 
(2) 

t  u\  a  f  F{v)dv 

v-r 

where 

LJ 

(=k 

and  up  is  the  plasma  frequency,  defined  by 


and  F(v)  ia  the  velocity  distribution  function  normalized  to 


There  are  in  general  four  solutions  u(k)  of  the  dispersion  relation,  given  real  k, 
(iwrj+wjj,)  and  (±uv,  These  were  found  numerically  using  Muller’s  method, 

and  one  of  each  pair  appears  in  Figure  1  for  the  example  V(  =  0.0.  For  small  enough 
k ,  one  pair  of  roots  has  a  positive  imaginary  part;  this  is  the  growth  rate  of  the 
mode  k.  Figure  2  shows  this  root  for  several  values  of  v*.  The  fastest  growth  rate  , 
the  mode  number  of  the  fastest  growing  mode  and  the  mode  number  at  which  the 
growth  rate  goes  through  zero  are  plotted  against  v(  in  Figure  3. 
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The  Penrooe  condition  gives  an  independent  determination  of  the  value  of  v( 
above  which  the  two  stream  distribution  is  stable  and  below  which  it  is  unstable. 
The  Penrose  condition  is  that 

for  stability  [4j,  where  /(f)  is  a  local  minimum  in  the  distribution  function,  /(v), 
normalized  to 

/  f(v)dv  =  N. 

For  the  two  stream  distribution  (Eq.(  2))  the  local  minimum  is  at  f  =  0.  A  numerical 
computation  using  Simpson’s  rule  showed  the  integral  in  Eq.(  5)  to  be  positive  for 
vt  >  0.77uo-  Thus  the  two  stream  distribution^r),  is  stable  if  vt  >  0.77«o,  and 
is  unstable  otherwise.  The  solutions  of  the  dispersion  relational;)  (Eq.(  3)), are 
consistent  with  this  result;a/(Ar)  shows  no  growing  modes  for  vt  >  0.77vo,  but  does 
show  growing  modes  for  v<  <  0.77«q. 

Simulation 

ESI  is  a  one  dimensional,  electrostatic  particle  simulation  code  developed  by  A. 
Bruce  Langdon(5].  Given  the  initial  positions  and  velocities  of  a  set  of  charged 
particles,  ESI  solves  for  the  electric  fields  on  grid  points  (in  this  case,  128)  and 
uses  the  fields  to  compute  new  particle  positions  and  velocities.  This  procedure  is 
repeated  each  time  step.  The  code  uses  periodic  boundary  conditions;  that  is,  if  a 
particle  leaves  one  end  of  the  system  it  will  reenter  with  the  same  velocity  at  the 
other  end.  The  length  of  the  system  used  is 

,  _  20*uo 
|  —  > 
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and  the  time  increment  is 


Input  parameters  instruct  ESI  to  load  the  particles  (typically  8192)  into  two  beams 
and  within  each  beam  to  distribute  the  velocities  in  a  Maxwellian  according  to  the 
thermal  velocity  specified. 

Linear  Results 

Logrithmic  plots  showing  the  growth  of  electrostatic  field  energy  in  the  various 
modes  are  generated  by  the  simulation.  Figure  4  contains  such  plots  for  rapidly 
growing  modes  for  several  values  of  vt.  The  theoretical  growth  rates  have  been 
sketched  in  for  comparison.  The  linear  nature  of  the  growth  is  apparent  in  the  cold 
cases,  and  agrees  with  the  theoretical  predictions. 

In  the  warmer  cases,  an  additional  effect  makes  such  a  comparison  more  diffi¬ 
cult.  There  seems  to  be  an  additional  contribution  to  the  growth  of  the  electrostatic 
energy  which  grows  nonlinear ly  and  with  a  much  greater  level  of  high  frequency  os¬ 
cillations  than  seen  in  the  cold  cases. The  effect  becomes  more  and  more  dominant  as 
the  value  of  v(  is  increased.  It  has  been  suggested  that  because  there  is  an  increasing 
difference  in  velocity  between  the  individual  particles  in  the  tail  of  the  Maxwellian 
as  the  temperature  is  increased,  perhaps  an  increasing  level  of  multibeaming  be¬ 
tween  these  tail  particles  is  being  observed  [6].  Two  experiments  were  conducted  to 
test  this  hypothesis.  The  number  of  particles  was  increased  in  order  to  reduce  the 
velocity  spread  between  particles  in  the  tail  of  the  distribution,  resulting  in  a  drop 
in  the  level  of  the  high  frequency,  nonlinear  contribution  to  the  electrostatic  energy. 
Another  test  was  to  observe  the  evolution  of  the  electrostatic  energy  of  a  single 
nondrifting  Maxwellian.  The  energy  did  grow  in  this  case,  and  the  growth  looked 
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very  similar  to  the  corresponding  (same  vt)  two  stream  case.  For  large  values  of 
V|,  the  multibeaming  instability  between  the  particles  in  the  tail  of  the  distribution 
may  make  a  larger  contribution  to  the  field  energy  than  the  two  stream  instability 
of  the  two  primary  beams.  ThiB  undesired  contribution  to  the  field  enery  can  be 
reduced  by  increasing  the  number  of  particles  in  the  simulation  (see  Figure  5),  but 
for  practical  numbers  of  particles  the  effect  still  obscures  the  linear  growth  for  some 
values  of  v(.  This  study  did  not  go  deeper  into  resolving  the  multibeaming  problem. 

Nonlinear  Results 

As  the  system  evolves  in  time,  depending  on  the  initial  temperature,  the  two  streams 
may  remain  distinct,  or  particle  mixing  may  take  place  until  the  original  beams  are 
no  longer  distinguishable.  In  the  cold  case  the  instability  is  severe,  and  the  two 
stream  character  is  lost.  Figure  6a  shows  plots  of  the  total  distribution  function, /(v), 
for  t>t  =s  0  at  intervals  of  twenty  time  steps.  Figure  6b  shows  plots  of  the  distribution 
of  velocities  of  the  particles  from  just  one  of  the  original  beams  at  the  same  time 
intervals.  After  the  seventh  plot  in  Figure  6b,  the  mean  velocity  has  decreased  to 
zero  from  its  initial  value  «o  =  1.  For  the  very  warm  example  in  Figure  7  (vt  =  0.7uo, 
still  unstable  by  the  Penrose  criterion)  the  instability  has  only  a  minor  effect  on  the 
distribution  and  the  two  stream  character  is  retained. 

In  some  of  the  intermediate  cases  an  additional  phenomenon  occurs.  The  results 
of  a  simulation  for  vt  =  0.1  are  shown  in  Figure  8.  The  total  distribution  oscillates 
between  a  double  peaked  and  single  peaked  curve  centered  at  v  =  0.  This  is  due 
to  the  spiralling  of  the  streams  in  phase  space  as  seen  in  the  sequence  of  Figure 
9.  Eventually,  the  instability  saturates,  oscillation  is  no  longer  observed  and  the 
distribution  function  /(v)  is  left  with  a  single  hump. 
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In  order  to  study  the  final  state  quantitatively,  a  single  parameter  is  introduced 
which  indicates  the  form  of  the  distribution  function.  Specifically,  No,  the  number 
of  particles  with  velocity  equal  to  zero  (in  practice,  the  number  of  particles  with 
velocity  in  the  range  -Au  <  v  <  Au,  where  Au  <  vq)  gives  a  measure  of  the  size 
of  the  dip  in  the  distribution  function  at  v  =  0  and  thus  a  measure  of  the  extent  to 
which  the  two  streams  have  remained  distinct.  Given  the  total  number  of  particles 
and  the  total  kinetic  energy  in  the  system,  a  single  Maxwellian  centered  at  t>  =  0  can 
be  completely  specified,  and  No  for  this  theoretical  limit  can  be  computed  from  the 
simulation  parameters,  and  is  called  Nom.  These  two  parameters  (total  number  of 
particles  and  kinetic  energy)  were  also  used  to  determine  a  two  stream  distribution 
with  peaks  at  v  =  ±1.  which  is  marginally  stable  by  the  Penrose  criterion,  allowing 
a  theoretical  estimate  of  No  for  marginal  stability,  called  N0p. 

If  the  instability  completely  randomized  the  particle  positions  in  phase  space, 
the  expected  final  distribution  of  velocities  would  be  a  single  Maxwellian  centered 
at  v  =  0.  Since  a  single  Maxwellian  has  no  dip  at  v  =  0.,  this  randomization  would 
result  in  an  upper  limit  on  the  value  of  No.  Alternatively,  if  the  instability  saturated 
at  the  point  of  marginal  stability,  a  minimum  value  of  No  would  result.  The  values 
of  No  found  in  the  simulations  have  been  compared  with  these  two  limits. 

Assuming  that  the  form  of  the  distribution  function  in  the  saturation  state  is 
bimaxwellian,  then  the  Penrose  condition  for  marginal  stability  is  vtp  =  0.77vo p, 
where  vtp  is  the  thermal  velocity  associated  with  the  two  beams.  Then  No  for  the 
Penrose  stable  case  (Nop)  is 


No,= 


2A  uN 
y/itK/m 


e~2  (rW)*, 


(6) 


and  for  the  single  Maxwellian, 


Nom  — 


2AuN 

y/irK/m 


(7) 


( K  is  the  average  kinetic  energy  per  particle,  and  N  is  the  number  of  particles  per 
beam.) 

Plots  of  Nq  vs.  time  for  various  initial  temperatures  are  shown  in  figure  10. 
The  values  of  No,  and  Nom  have  been  sketched  in  for  comparison.  As  expected,  No 
decreases  with  respect  to  the  theoretical  values  as  the  beams  are  warmed.  In  the 
colder  cases,  it  is  somewhere  between  the  two  extremes,  and  in  the  warmer  cases  it 
is  near  No,-  The  large  oscillations  seen  in  No,  particularly  for  v(  =  0.1,  are  due  to 
the  spiralling  of  the  streams  in  phase  space  as  described  previously. 


Conclusion 


The  dispersion  relation  for  the  instability  changed  smoothly  as  the  thermal  velocity 
was  increased.  The  growth  rates  from  the  simulations  agreed  with  the  theoretical 
values,  to  the  extent  to  which  they  could  be  obtained.  The  final  shape  of  /(v)  varied 
gradually  with  v(  from  a  single  peaked  distribution  for  the  cold  case  to  an  almost 
unchanged  double  peaked  distribution  in  the  warm  case.  No  sudden  transitions  were 
observed;  i.  e.,  there  appears  to  be  no  special  value  of  vt/vo- 
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APPENDIX  A-  Calculation  of  iVom>  the  number  of  particles  with  velocities 
in  the  range  -Au  <  v  <  Au  for  a  single  Maxwellian  centered  on  v  =  0  with  the 
same  total  number  of  particles  and  total  kinetic  energy  as  the  initial  two  stream 
distribution  used  in  the  particle  simulations. 


In  order  to  compute  Nom » it  is  necessary  to  compute  the  root  mean  square  velocity, 
vm  of  the  particles  in  the  single  Maxwellian.  vm  can  be  found  by  equating  the  total 
energy  of  the  particles  in  the  single  Maxwellian  to  the  total  particle  energy  in  the 
two  stream  distribution: 

i(2iV)mv^  =  i(2jV)m(t>g  +  vf).  (8) 

Thus, 

vl  =  vo  +  vt-  (9) 


The  normalization  for  the  distribution  is  found  by  integrating  over  velocities  to 
get  the  total  number  of  particles: 

jrl*  =  2*-  (l0) 

N  is  the  number  of  particles  in  each  beam.  The  normalization  constant  is: 


By  integrating  the  normalized  distribution  function  between  -Au  and  Au,  the 
number  of  particles  with  velocity  between  these  limits  can  be  computed: 


dv. 


(12) 
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Since  Au  <  vq,  the  integral  may  be  approximated  by  multiplying  the  length  of  the 
interval, 2Au,  by  the  value  of  the  integrand  evaluated  at  v  =  0.,  resulting  in  the  final 


expression 

i V0m  a  (13) 

Vm  V  7T 

The  average  kinetic  energy  per  particle  is  thus 

K  =  -mvL. 

2  m 

APPENDIX  B-  Calculation  of  IVop,  the  number  of  particles  with  velocities 
in  the  range  — Au  <  v  <  Au  for  a  bimaxwellian  centered  on  v  =  0  with  vt  =  0.77 
and  with  the  same  total  number  of  particles  and  total  kinetic  energy  as  the  initial 
two  stream  distribution  used  in  the  particle  simulations. 


This  calculation  is  analogous  to  that  for  Nom-  The  total  energy  of  the  particles  in  a 
marginally  Penrose  stable  bimaxwellian  is 


NK=\m{2N)  (v?p  +  v*p) 

where  K  is  again  the  average  kinetic  energy  per  particle,  t/<>p  is  the  beam  velocity  and 
vtp  is  the  thermal  velocity.  For  the  marginally  stable  case,  vlp  =  0.77vop-  Equating 
this  energy  to  the  total  energy  given  of  the  particles  in  the  simulations, 


Thus, 


(21V)  (21V)  (v%  +  v?)  . 

2  _  v0  +  v?  _  *>0  +  vt 

*  (l.  +  (0.77)2)  1593  ’ 
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\  a  15.9 


t s  35.9 


Figure  5:  Field  energy  vt.  time  for  vt  =  0.5  for  the  cues  (a)  16384  particles  per  beam  and  (b) 
65536  particles  per  beam. 


are  3:  Dispersion  relation  numerical  solutions:  plots  of  (a)  maximum  growth  rate  vs.  thermal 
>city  of  beams,  (b)  wave  number  of  fastest  growing  mode  vs.  thermal  velocity,  and  (c)  the  mode 
iber  at  which  the  growth  rate  goes  through  zero  vs.  the  thermal  velocity.  is  normalized  to 
k  is  normalized  to  vq/up  and  vt  is  normalized  to  vq. 
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and  (f),  and  the  positive  root  is  shown  in  (b)  and  (c).  Note  that  both  positive  and  negative  roots 
exist  for  all  cases,  but  that  only  one  root  was  plotted. 
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The  appropriate  Penrose  stable  distribution  is  thus 


Normalizing  to  find  the  constant  B, 

f°°  fp(u)du  =  2 N 
J— 00 

gives 


B  = 


N 


vtpV2*‘ 

Integrate  fp(u)  between  — Au  and  Au  to  find  Nop-. 


(**  [  2*?,  J+“P[  2^  J  j  ’ 


or 


2Au  N 

Nop  SJ - 7==  exp 

v»fV2x 


"bp 

2v?p 


=  0.43 


2AuN 

y/xK/m' 
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same  initial  parameters  is  also  indicated. 


